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Abstract 

Evolutionary graph theory studies the evolutionary dynamics of populations 
structured on graphs. A central problem is determining the probability that 
a small number of mutants overtake a population. Currently, Monte Carlo 
simulations are used for estimating such fixation probabilities on general di- 
rected graphs, since no good analytical methods exist. In this paper, we 
introduce a novel deterministic framework for computing fixation probabil- 
ities for strongly connected, directed, weighted evolutionary graphs under 
neutral drift. We show how this framework can also be used to calculate 
the expected number of mutants at a given time step (even if we relax the 
assumption that the graph is strongly connected), how it can extend to other 
related models (e.g. voter model), how our framework can provide non-trivial 
bounds for fixation probability in the case of an advantageous mutant, and 
how it can be used to find a non-trivial lower bound on the mean time to 
fixation. We provide various experimental results determining fixation prob- 
abilities and expected number of mutants on different graphs. Among these, 
we show that our method consistently outperforms Monte Carlo simulations 
in speed by several orders of magnitude. Finally we show how our approach 
can provide insight into synaptic competition in neurology. 
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1. Introduction 

Evolutionary graph theory (EGT), introduced by [Hj, studies the prob- 
lems related to population dynamics when the underlying structure of the 
population is represented as a directed, weighted graph. This model has 
been applied to problems in evolutionary biology [28], physics [25], game 
theory [20], neurology [26.J, and distributes systems [13] . A central problem 
in this research area is computing the fixation probability - the probability 
that a certain subset of mutants overtakes the population. Although good 
analytical approximations are available for the undirected/unweighted case 
[U E], these break down for directed, weighted graphs as shown by [16J. As 
a result, most work dealing with evolutionary graphs rely on Monte Carlo 
simulations to approximate the fixation probability [221 II]- in this paper 
we develop a novel deterministic framework to compute fixation probability 
in the case of neutral drift (when mutants and residents have equal fitness) 
in directed, weighted evolutionary graphs based on the convergence of "ver- 
tex probabilities" to the fixation probability as time approaches infinity. We 
then show how this framework can be used to calculate the expected number 
of mutants at a given time, how the framework can be modified to do the 
same for related models, how it can provide non-trivial bounds for fixation 
probability in the case of an advantageous mutant, and how it can provide 
a non-trivial lower bound on the mean time to fixation. We also provide 
various experiments that show how our method can outperform Monte Carlo 
simulations by several orders of magnitude. Additionally, we show that the 
results of this paper can provide direct insight into the problem of synaptic 
competition in neurology. 

Our method also fills a few holes in the literature. First, it allows for 
deterministic computation of fixation probability when there is an initial 
set of mutants - not just a singleton (the majority of current research on 
evolutionary graph theory only considers singletons). Second, it allows us 
to study how the mutant population chang function of time. Third, 

we show (by way of rigorous proof) that fixation probability, under the case 
of neutral drift is a lower bound for the case of the advantageous mutant - 
confirming simulation observations by [15]. Fourth, we show (also by way of 
rigorous proof) that fixation probability under neutral drift is additive (even 
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for weighted, directed graphs), which extends the work of £6] that proved 
this for undirected/unweighted graphs. Fifth, we provide a non-trivial lower 
bound for the computation of mean time to fixation in the general case - 
which has only previously been explored for well-mixed populations [2] and 
special cases of graphs [3]. 

This paper is organized as follows. In Section [2] we review the original 
model of Lieberman et al., introduce the idea of "vertex probabilities" and 
show how they can be used to find the fixation probability. We then show how 
this can be used to determine the expected number of mutants at a given 
time in Section [3} This is followed by a discussion of how the framework 
can be extended to other update rules in Section [4] and then for bounding 
fixation probability in the case of an advantageous mutant in Section |5j 
We then discuss how our approach can be adopted to bound mean time to 
fixation in Section |6j We use the results of the previous sections to create 
an algorithm for computing fixation probability and introduce a heuristic 
technique to significantly decrease the run-time. The algorithm and several 
experimental evaluations are described in Section [7| In Section [8j we show 
how our framework can be applied to neurology to gain insights into synaptic 
competition. Finally, we discuss related work in Section [9] and conclude. 

2. Directly Calculating Fixation Probability 

The classic evolutionary model known as the Moran Process is a stochas- 
tic process used to describe evolution in a well- mixed population |18j . All 
the individuals in the population are either mutants or residents. The aim of 
such work was to determine if a set of mutants could take over a population 
of residents (achieving "fixation"). In [14], evolutionary graph theory (EGT) 
is introduced, which generalizes the model of the Moran Process by speci- 
fying relationships between the N individuals of the population in the form 
of a directed, weighted graph. Here, the graph will be specified in the usual 
way as G = (V, E) where V is a set of nodes (individuals) and E C V x V. 
In most literature on evolutionary graph theory, the evolutionary graph is 
assumed to be strongly connected. We make the same assumption and state 
when it can be relaxed. 

For any node i, the numbers k\^, k„l t are the in- and out- degrees respec- 
tively. We will use the symbol N to denote the sized of V. Additionally, 
we will specify weights on the edges in set E using a square matrix denoted 
W = [wij] whose side is of size N. Intuitively, Wij is the probability that 
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member of the population j is replaced by % given that member % is selected. 
We require ^ ■ Wij = 1 and that G E iff u>jj > 0. If for all i, j, we have 

Wij = l/kol t , then the graph is said to be "unweighted." If for all G E, 
we have (J,i) G the graph is said to be "undirected." Though our results 
primarily focus on the general case, we will often refer to the special case 
of undirected/unweighted graphs as this special case is quite common in the 
literature jH |6]. 

In this paper we will often consider the outcome of the evolutionary pro- 
cess when there is a set of initial mutants as opposed to a singleton. Hence, 
we say some set (often denoted C) is a configuration if that set specifies the 
set of mutants in the population (all other members in the population then 
are residents). We assume all members in the population are either mutants 
or residents and have a fitness specified by a parameter r > 0. Mutants 
have a fitness r and residents have a fitness of 1. At each time step, some 
individual i G V is selected for "birth" with a probability proportional to its 
fitness. Then, an outgoing neighbor j of % is selected with probability 
and replaced by a clone of i. Note if r = 1, we say we are in the special case 
of neutral drift. 

We will use the notation Pv>,t to refer to the probability of being in 
configuration V after t timesteps and Py,t\c to be the probability of being in 
configuration V at time t conditioned upon initial configuration C. Perhaps 
the most widely studied problem in evolutionary graph theory is to determine 
the fixation probability. Given set of mutants C at time 0, the fixation 
probabilty is defined as follows. 

F c = lim t ^Py AC (1) 

This is the probability that an initial set C of mutants takes over the en- 
tire population as time approaches infinity. Similarly, we will use the term 
the extinction probability, Fc, to be lim t -^ooP(D t t\c- If the graph if strongly 
connected, then we have Fc + Fc = 1. Hence, for a strongly connected 
graph, a mutant either fixates or becomes extinct. Typically, this problem 
is studied using Monte Carlo simulation. This work uses the idea of a ver- 
tex probabilities to create an alternative to such an approach. The vertex 
probability is the probability that a certain vertex is a mutant at a certain 
time given an initial configuration. For vertex % at time t, we denote this as 
Pi,t\c- Often, for ease of notation, we shall assume that the probabilities are 
conditioned on some initial configuration and drop the condition, writing P i t 
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instead of Pi,t\c- We note that P iyt can be expressed in terms of probabilities 
of configurations as follows. 

P itt = Yl p vi (2) 

V'£2 V 

s.t. ieV 

Viewing the probability that a specific vertex is a mutant at a given time 
has not, to our knowledge, been studied before with respect to evolutionary 
graph theory (or in related processes such as the voter model). The key 
insight of this paper is that studying these probabilities sheds new light on 
the problem of calculating fixation probabilities in addition to providing other 
insights into EGT. For example, it is easy to show the following relationship. 

Proposition 1. Let V be a subset of V and t be an arbitrary time point. 
Iff for all i G V , P itt = 1 and for all i V , P iti = 0, then Pv,t = 1 and for 

all V" G 2 V S.t. V" V, Pyn-t = 0. 

It is easy to verify that Fc > iff Wi G V, Umt->ooPi,t > 0. Hence, in this 
paper, we shall generally assume that lim t ->ooPi,t > holds for all vertices i 
and specifically state when it does not. As an aside, for a given graph, this 
assumption can be easily checked: simply ensure for j G V — C that exists 
some ieC s.t. there is a directed path from i to j. 

Now that we have introduced the model and the idea of vertex proba- 
bilities we will show how to leverage this information to compute fixation 
probability. It is easy to show that as time approaches infinity, the vertex 
probabilities for all vertices converge to the fixation probability when the 
graph if strongly connected. 

Theorem 1. Vi, lim t -yoo P%,t\c = Fc 

Now let us consider how to calculate for some i and t. For t — 0, 
where we know that we are in the state where only vertices in a given set are 
mutants, we need only appeal to Proposition[T]- which tells us that we assign 
a probability of 1 to all elements in that set and otherwise. For subsequent 
timesteps, we have developed Theorem [2] shown next (the proof of which is 
included in the supplement). 
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Theorem 2. 

Pi,t = Pi,t-1 + w i i fe- 1 ' %*)[(?.*-!) ~ ' %t)|(i,t-l)) 

(^Cj'i*) I ^ s ^ e probability that j is picked for reporduction at time t given 
that % was a mutant at time t — 1.) 

We believe that a concise, tractable analytical solution for S(j t t)\u,t-i) is un- 
likely. However, for neutral drift (r = 1), these conditional probabilities are 
trivial - specifically, we have for all i, j, t, S(j ; t)\(i,t-X) — as this probability 
of selection is independent of the current set of mutants or residents in the 
graph. Hence, in the case of neutral drift, we have the following: 

Pi,t = Pi,t-t + £ f ■ (P 3 ,t-i - P,t-i) (3) 

Studying evolutionary graph theory under neutral drift was a central theme 
in several papers on EGT in the past few years [HI US] as it provides an 
intuition on the effects of network topology on mutant spread. In Section [5] 
we examine the case of the advantageous mutant (r > 1). Neutral drift allows 
us to strengthen the statement of Equation [T] to a necessary and sufficient 
condition - showing that when the probabilities of all nodes are equal, then 
we can determine the fixation probability. 

Theorem 3. Assuming neutral drift (r = 1), given initial configuration C 
with fixation probability Fq, if at time t the quantities P%,t\c are equal (for all 
i G V ), then they also equal Fc- 

Therefore, under neutral drift, we can determine fixation probability when 
Equation [3] causes all P^t's to be equal. We can also use Equation [3] to find 
bounds on the fixation probability for some time t by the following result 
that holds for any time t under neutral drift. 

min P i}t <F C < max P i t (4) 

i i 

Under neutral drift, we can show that fixation probability is additive for 
disjoint sets. Broom et al. proved a similar result the a special case of undi- 
rected/unweighted evolutionary graphs [6j. However, our proof (contained in 
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the supplement) differs from theirs in that we leverage Equation [3] Further, 
unlike the result of Broom et al., our result applies to the more general case 
of weighted, directed graphs. 

Theorem 4. When r = 1 for disjoint sets C, D C V, F c + F D = F C ud- 



3. Calculating the Expected Number of Mutants 

In addition to allowing for the calculation of fixation probability, our 
framework can also be used to observe how the expected number of mutants 
changes over time. We will use the notation Ex^ to denote the expected 
number of mutants at time t given initial set C. Formally, this is defined 
below. 

Exg' = $>,« (5) 

i£V 

Unlike fixation probability, which only considers the probability that mu- 
tants overtake a population, Ex|J' provides a probabilistic average of the 
number of mutants in the population under a finite time horizon. For exam- 
ple, is has been noted that graph structures which amplify fixation normally 
also increase time to absorption 0, |2T] . Hence, finding the expected number 
of mutants may be a more viable topic in some areas of research where time 
is known to be limited. Following from Equation [3] where we showed how to 
compute Pi tt for each node at a given time, we have the following relationship 
concerning the expected number of mutants at a given time under neutral 
drift. 



Ex? = Exg- 1} + ^ - i £ E P ^ (8) 

i&V (j,i)eE 

Based on Equation [6j we notice that for r = 1, at each time-step, the 
number of expected mutants increases by at most the average fixation prob- 
ability and decreases by a quantity related to the average "temperature." 
The temperature of vertex i (denoted Tj) is defined for a given node is the 
sum of the incoming edge weights p3J: T% = Yl,j w ji- Intuitively, nodes with 
a higher temperature change more often between being a mutant and being a 
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resident than those with lower temperature. Re- writing Equation [6] in terms 
of temperature we have the following: 

Ex» = Exg-n^ l^.P^ (7) 

iev 

Hence, if the preponderance of high temperature nodes are likely to be 
mutants, then most likely the average number of mutants will decrease at the 
next time step. We also note that Theorem |2j Equation [3j and Equation [6] 
do not depend on the assumption that the underlying graph is strongly con- 
nected. Therefore, as such is the case, we can study the relationship of time 
vs. expected number of mutants for any evolutionary graph (under neutral 
drift). This could be of particular interest to non-strongly connected evolu- 
tionary graphs that may have trivial fixation probabilities (i.e. 1 or 0) but 
may have varying levels of mutants before achieving an absorbing state. 

4. Applying the Framework to Other Update Rules 

The results of the last two sections not only apply to the original model of 
|14j . but several other related models in the literature. Viewing an evolution- 
ary graph problem as a stochastic process, where the states represent different 
mutant-resident configurations, it is apparent that the original model spec- 
ifies the transition probabilities. However, there are other ways to specify 
the transition probabilities known as update rules. Several works address 
different update rules pfl |25l [15]. Overall, we have identified three major 
families of update rules - birth-death (a.k.a. the invasion process) where the 
node to reproduce is chosen first, death-birth (a.k.a. the voter model) where 
the node to die is chosen first, and link dynamics, where an edge is chosen. 
We summarize these in Table [U 

We have already shown how our methods can deal with the original model 
of Lieberman et al., often referred to as the Birth-Death (BD) process. In 
this section, we apply our methods to the neutral-drift (non-biased) cases of 
death-birth and link-dynamics. In these models, the weights of the edges is 
typically not considered. Hence, in order to align this work with the majority 
of literature on those models, we will express vertex probabilities in terms of 
node in-degree (&•« ) an d the set of directed edges (E). We note that these 
results can be easily extended to a more general case with an edge-weight 
matrix as we used for the original model of EGT. 
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Table 1: Different families of update rules. 



Update Rule 


Intuition 


Birth-Death (BD) 

(a.k.a. Invasion Process (IP)) 


(1) Node i selected 

(2) neighbor of i, node j selected 

(3) Offspring of i replaces j 


Death-Birth (DB) 

(a.k.a. Voter Model (VM)) 


(1) Node i selected 

(2) neighbor of i, node j selected 

(3) Offspring of j replaces i 


Link Dynamics (LD) 


(1) Edge selected 

(2) The offspring of one node in the 
edge replaces the other node 



4-1. Death- Birth Updating 

Under the death birth model (DB), at each time step, a vertex i is selected 
for death. With a death-bias (DB-D), it is selected proportional to the inverse 
of its fitness, with a birth-bias (DB-B) it is selected with a probability 1/N, 
which is also the probability under neutral drift. Then, an incoming neighbor 
(j) is selected either proportional to the fitness of all incoming neighbors 
(birth-bias), or with a uniform probability (in the case death-bias or neutral 
drift). The selected neighbor then replaces i. Here, we compute P^t under 
this dynamic with r — 1. 

Pi,t = (1 - N' 1 )^ + (Nk^y 1 Yl Pj*-i ( 8 ) 

We note that the proof of convergence still holds for death-birth - that 
is for some time t, Vi, the value P^ t is the same, then P^ t = Fq. Further, 
Theorem [4] holds for DB under neutral drift as well, specifically, for disjoint 
sets C, D C V, F c + P D = Pcud- 

4-2. Link- Dynamics 

With link dynamics (LD), at each time step an edge is selected 

either proportional to the fitness of i or the inverse of the fitness of j. It 
has previously been shown that LD under birth bias is an equivalent process 
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to LD with a death bias [IB] . Under neutral drift, the probability of edge 
selection is 1/\E\ (where \E\ is the cardinality of set E). Then, i replaces j. 
Now, we compute P^t under this dynamic with r = 1. 

= (i - A^iEr 1 )^-! + ]4 E p ^ ( 9 ) 

Again, convergence and additivity of the fixation probability still hold 
under link dynamics just as with BD and DB. 

5. Bounding Fixation Probability for r > 1 

So far we have shown how our method can be used to find fixation proba- 
bilities under the case of neutral drift. Here, we show how our framework can 
be useful in the case of an advantageous mutant (when the value for r, the 
relative fitness, is greater than 1). First, we show that our method provides 
a lower bound. We then provide an upper bound on the fixation probabil- 
ity that can be used in conjunction with our framework when studying the 
case of the advantageous mutant. We note that certain parts of these proofs 
are specific for diffent update rules, and we identify them using the abbre- 
viations from the last section (DB-D, DB-B, and LD). The update of the 
original model of [2] is known as the "birth-death" model and abbreviated 
BD. If the fitness bias is on a birth event, we denote it as BD-B and if the 
bias is on a death event we denote it as BD-D. 

Naoki Masuda observes experimentally (through simulation) that the fix- 
ation probability computed with neutral drift appears to be a lower bound 
on the fixation probability for an advantageous mutant [15J. We were able 
to prove this result analytically - the proof is included in the supplementary 
materials. 

Theorem 5. For a given set C , let be the fixation probability under 

(r) 

neutral drift and Fq be the fixation probability calculated using a mutant 
fitness r > 1. Then, under BD-B, BD-D, DB-B, DB-D, or LD dynamics, 

r c — r c ■ 

This proof leads to the conjecture that r' > r implies Fq ^ > Fq . How- 
ever, we suspect that proving this monotonicity property will require a dif- 
ferent technique than used in Theorem [5] Next, to find an upper bound 
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that corresponds with the lower bound above, we use the proof technique in- 
troduced in [9], to obtain the following non-trivial upper bounds of fixation 
probability for individual nodes in various update rules. 



BD-B: F {i} < r(r + J2^ji)~ 



(10) 
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BD-D: F {i} < (y ^ ) 



(11) 



DB - B : F {i} < y rwij{\ - + rwij) 



-i 



(12) 



j 



DB-D: F {i} < ry Wij 



(13) 
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6. A Lower Bound for Mean Time to Fixation 

Another important, although less-studied problem with respect to evolu- 
tionary graph theory is the mean time to fixation - the average time it takes 
for a mutant to take over the population. Closely related to this problem 
are mean time to extinction (average time for the resident to take over) and 
mean time to absorption (average time for either mutant or resident to take 
over). This has been previously studied under the original Moran process 
for well mixed populations |2j as well as some special cases of graphs [3]. 
However, to our knowledge, a general method to compute these quantities 
(without resorting to the use of simulation) have not been previously studied. 
Here we take a "first step" toward developing such a method by showing how 
the techniques introduced in this paper can be used to compute a non-trivial 
lower bound for mean time to fixation (and easily modified to bound mean 
time to extinction and absorption). 

Let F t \c be the probability of fixation at time t. Therefore, F t \c — F t -i\c 
is the probability of entering fixation at time t. The symbol tc is the mean 
time. By the results of [2], we have the following: 

Theorem 6. 



-j^y*- (Ft\c - Ft-i\c) 
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Our key intuition is noticing that at each time step t, Fqc < min^ P i t . 
From this, we use the accounting method to provide a rigorous proof for the 
following theorem that provides a non-trivial lower-bound for the mean time 
to fixation. This result can be easily modified for mean time to extinction 
and absorption as well. 

Theorem 7. ^ £ t =i t- (P min ,*-iW-i) < ^ E t =i *■ (F t{c - F t _ 1{c ) Where 
Pmm,t — minj Pi t . 

7. Algorithm and Experimental Evaluation 

We leverage the finding of the previous sections in Algorithm [TJ As described 
earlier, our method has found the exact fixation probability when all the 
probabilities in Uj{Pi,*} (represented in the pseudo-code as the vector p) are 
equal. We use Equation [4] to provide a convergence criteria based on value 
e, which we can prove to be the tolerance for the fixation probability. 

Proposition 2. Algorithm^ returns the fixation probability Fq within ±e. 

Our novel method for computing fixation probabilities on strongly con- 
nected directed graphs allows us to compute near-exact fixation probabilities 
within a desired tolerance. The running time of the algorithm is highly de- 
pendent on how fast the vertex probabilities converge. In this section we 
experimentally evaluate how the vertex probabilities in our algorithms con- 
verge. We also provide results from comparison experiments to support the 
claim that Algorithm^ AC C finds adequate fixation probabilities order of 
magnitudes faster than Monte Carlo simulations. We also show how the al- 
gorithm can be used to study the expected number of mutants as well as 
bound mean time to fixation. 

7.1. Convergence of Vertex Probabilities 

We ran our algorithm to compute fixation probabilities on randomly 
weighted and strongly connected directed graphs in order to experimentally 
evaluate the convergence of the vertex probabilities. We generated the graphs 
to be scale-free using the standard preferential attachment growth model [3] 
and randomly assigned an initial mutant node. We replaced all edges in the 
graph given by the growth model with two directed edges and then randomly 
assigned weights to all the edges. 
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Algorithm 1 - Our Novel Solution Method to Compute Fixation Probabil- 
ities 

Input: Evolutionary Graph (N, V, W), configuration C C V, natural number 

R > 0, and real number e > 0. 

Output: Estimate of fixation probability of mutant. 

1: pi is the ith position in vector p corresponding with vertex i G V. 
2: Set Pi = 1 if i G C and Pi = otherwise. {As per Proposition □]} 
3: q <r- p {q will be p from the previous time step.} 

4: T <r- 1 

5: while r > e do 

6: for i G V {This loop carries out the calculation as per Equation ^ 
do 

7: sum ^— 

8: m <- {j G V\wji > 0} 

9: for j G m do 
10: sum = sum + Wji ■ (qj — q^) 

11: end for 
12: pi qj + 1/N ■ sum 
13: end for 
14: q ^— p 

15: r <— (1/2) • (maxp — minp) {Ensures error bound based on Equation [4]f 

16: end while 

17: return (minp) + r 
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Figure 1: Left: Convergence of the minimum (MinP), maximum (MaxP), and average 
(AvgP) of vertex probabilities towards the final fixation probability as a function of our 
algorithm's iterations t for a graph of 100 nodes. Right: Average speedup (on a log scale) 
for finding fixation probabilities achieved by our algorithm vs Monte Carlo simulation for 
graphs of different sizes. 



To compare Algorithm [T] with the Monte Carlo approach, we should 
set the parameter R in that algorithm to be comparable with e in Algo- 
rithm [TJ As e is the provable error of a solution to Algorithm [TJ Based on 
the commonly-accepted definition of estimated standard error from statis- 
tics, we can obtain the estimated standard error for the solution returned by 
Monte Carlo approach with the following expression (where R is the number 
of simulation runs). 



MLzjg) (14) 

R-l y 1 



We can use Equation 14 to estimate the parameter R for the Monte 



Carlo approach as follows. We set e equal to the estimated standard error 



as per Expression [14] and manipulate it algebraically. This gives us R ~ 
S ^ 1 ^ + 1 where S is the solution to Algorithm [lj e is the input parameter 
for Algorithm [T] and R is the number of simulation runs in the Monte Carlo 
approach that we estimate to provide a comparable error bound. We also 
note, that as the vertex probabilities converge, the standard deviation of the 
p vector in Algorithm [T] could be a potentially faster convergence criteria. 
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Note that using standard deviation of p and returning the average vertex 
probability would no longer provide us of the guarantee in Proposition |2j 
however it may provide good results in practice. The modifications to the 



algorithm would be as follows: line 15 would be r f- st.dev(p) and line 17 



would be return avg(p). We will refer to this as Algorithm^ with alternate 
convergence criteria or Algorithm^ AC C for short. 

Figure [l] (left) shows the convergence of the minimum, maximum, and the 
average of vertex probabilities towards the final fixation probability value for 
a small graph of 100 nodes. We can observe that the average converges to 
the final value at a logarithmic rate and much faster than the minimum and 
maximum vertex probability values. This suggests that while Algorithm^ 
ACC does not give the same theoretical guarantees as Algorithm^ it is much 
preferable for speed since the minimum and maximum vertex probabilities 
take much longer to converge to the final solution than the average. The 
fact that the average of the vertex probabilities is much preferable as a fast 
estimation of fixation probabilities is supported by the logarithmic decrease 
of the standard deviation of vertex probabilities (see Figure [2]). Convergences 
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for other and larger graphs are not shown here but are qualitatively similar 
to the relative convergences shown in the provided graphs. 



7.2. Speed Comparison to Monte Carlo Simulation 

In order to compare our method's speed compared to the standard Monte 
Carlo simulation method, we must determine how many iterations our algo- 
rithm must be run to find a fixation probability estimate comparable to that 
of the Monte Carlo approach. Thankfully, as we have seen, we can get a 
standard error on the fixation probability returned by the Monte Carlo ap- 



proach as per Equation 14 While we did not theoretically prove anything 
about how smoothly fixation probabilities from our methods approach the 
final solution, the convergences of the average and standard deviation as 
shown above strongly suggest that estimates from our method approach the 
final solution quite gracefully. In fact, in the following experiments, once 
our method has arrived at a fixation probability estimate within the stan- 
dard error of simulations, the estimate never again fell outside the window of 
standard error (although the estimate did not always approach the final es- 
timate monotonically) . This is in stark contrast to Monte Carlo simulations, 
from which estimations can vary greatly before the method has completed 
enough single runs to achieve a good probability estimate. 

We generated a number of randomly weighted and strongly connected 
directed graphs of various sizes on which we compare our solution method 
to Monte Carlo approximation of fixation probabilities. The graphs were 
generated as in our convergence experiments. For each graph of a different 
size, we generated a number of different initial mutant configurations. We 
found fixation probabilities both using Monte Carlo estimation with 2000 
simulation runs and our direct solution method, terminating when we have 
reached within the standard error of the Monte Carlo estimation. Since the 
average vertex probability proved to be such a good fast estimate of the true 
fixation probability, we used Algorithm^ AC C. 

Figure [I] (right) shows the speedup our solution provides over Monte 
Carlo simulation. Here speedup is defined as the ratio of the time it takes 
for simulations to complete over the time it takes our algorithm to find a 
fixation probability within the standard deviation. The often extremely low 
number of iterations needed by our algorithm to find fixation probabilities 
within the standard error of simulations may prompt the concern that the 
probabilities fall within this window so soon by mere chance. However, our 
experiments have shown that the fixation probability estimation found by 
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our algorithm at each iteration approaches the final fixation probability after 
termination smooothly at a logarithmic rate, asymptotically approaching the 
true fixation probability. While in this case the fixation probability estimate 
slightly crosses over the true fixation probability and then slowly approaches 
it again, none of the fixation probability estimates from our algorithm exited 
the window of standard error (from simulations) once they entered it. 

We can observe from our experiments that computing fixation proba- 
bilities using Monte Carlo simulations showed to be a very time-expensive 
process, highlighting the need for faster solution methods as the one we have 
presented. Especially for larger graph sizes, the time complexity of our so- 
lution to achieve similar results to Monte Carlo simulation has shown to be 
orders of magnitude smaller than the standard method. 

7. 3. Monitoring the Expected Number of Mutants 

As observed in Section [3j our method not only allows for the calculation 
of the fixation probability of a mutant, but also allows us to study how the 
expected number of mutants change over time. In this section, we present 
experimental results exploring the trajectory of the expected number of mu- 
tants over time on various undirected/unweighted graphs and under different 
initial mutant placement conditions. 

First, we note that the expected number of mutants (as time approaches 
infinity) in an unweighted/undirected graph with respect to a single initially 
infected vertex i can be computed by modifying the result of ||6] (for BD 
updating) to obtain the following. 

Urn Ex$ = - - Tx (15) 

Where is the average inverse of the degree for the graph. Hence, we 

can determine whether a node amplifies or suppresses selection by observing 
if lim^oo Exj*| is greater or less than 1 respectively: if ki < -^Ixy selection is 
amplified and if ki > ^pyr it is suppressed. We have used our algorithm to 
compare the trajectory of the expected number of mutants over time when 
the initial mutant is placed on amplifiers vs. suppressors under different 
graph topologies and BD updating. We note that similar comparisons can 
be obtained with our algorithm for the other update rules. We also note that 
by Theorem [5j an amplifier for BD (with no bias) will also be an amplifier 
for the (biased) BD-B and BD-D where r > 1. 
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Figure 3: Expected number of mutants over time starting with a single mutant placed 
on a graph for Barabasi- Albert preferential attachment (BAR), Erdos-Renyi (ERD), and 
Newmann-Watts-Strogatz small world (NWS) graphs. Lines are averages over 50 random 
graphs of each type. In the left graph, mutants are placed at the highest degree nodes, 
which are suppressors. In the right graph, mutants are placed at lowest degree nodes, 
which are amplifiers. 



Figure [3] shows the trajectories of the expected number of mutants over 
time on random [3] preferential attachment (BAR), pi)] (ER), and [19] small 
world graphs (NWS), each for when the initial mutant is placed on a sup- 
pressor (highest degree node of graph) and amplifier (lowest degree node of 
graph). Graphs are all of equal size at 100 nodes. We note that the highest 
degree nodes are especially strong suppressors on BAR graphs, less so for 
NWS graphs, and even less so for ER graphs. This makes sense when one 
considers the degree distribution of the different graph topologies, which are 
scale-free or power-law (P(k) ~ k~ 3 ) for BR, roughly Poisson-shaped for 
NWS, and relatively uniform for ER graphs. For lowest degree amplifiers, 
the expected number of mutants grows faster early on in Barabasi-Albert 
graphs, but it plateaus earlier than and is eventually surpassed by the slower 
growing expected number of mutants in the Erdos-Renyi, and Newmann- 
Watts-Strogatz graphs. Such insights into the evolutionary process may be 
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Figure 4: Expected number of mutants over time for an Erdos-Renyi graph of 100 nodes, 
with an extra muntant node (mn) and resident node(rn) with directed edges to conjoin 
and conjrn respectively. The value that the expected number of mutants converges to 
depends on the relative degrees of conjmn and conjrn, as shown in the legend. 



crucial in applications, e.g. when one may be more interested in achieving 
highest number of mutants in a short amount of time rather than highest 
number of mutants as t — > oo or vice versa. 

Finally, thus far we have only considered strongly connected graphs in 
which the vertex probabilities converge as t — > oo, but this is not the case 
for some non-strongly connected graphs. We have thus also investigated the 
expected number of mutants over time for some simple cases of such graphs. 
Consider a random graph that is strongly connected, and then have a resident 
node (rn) and mutant node (mn) connected with only directed edges into the 
strongly connected graph. Clearly, the vertex probabilities cannot converge, 
since Vt, P mn ,t = 1 and P rri)t = 0. Our experimental results in Figure [4] show 
however that while the vertex probabilities do not converge, the value for 
the expected number of mutants given by our algorithm seems to converge. 
What value the expected number of mutants converges to depends on the 
relative degrees of the nodes that the mutant node mn and resident node 
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rn connect to. We shall call these nodes con.mn and corurn, respectively. 
If k conmn ~ k conrn , the expected value of mutants converges at around 50% 
of the graph's nodes. If k conmn > k conrn , the expected value of mutants 
is less than 50% of the graph's nodes, and conversely, if k con mn < k con rn , 
it is greater. These results are intuitive because lower degree nodes are 
better spreaders under BD updating. These results are also interesting be- 
cause the expected value converges - even though the graphs are not strongly 
connected. By an examination of Equation |6j this convergence is possible. 
However, we have not proven that convergence always occurs. An interest- 
ing direction for future work is to identify under what conditions will the 
expected number of mutants converges in a non-strongly connected graph. 

7.4- Experimentally Computing the Lower Bound of the Mean Time to Fix- 
ation 

We also performed experiments to examine the lower bound on mean 
time to fixation (discussed in Section [6]) as compared to the average fixation 
time determined from simulation run. In doing so, we were able to confirm 
the lower-bound experimentally. We were able to use Algorithm [IJ-ACC to 
compute the lower bound with a few changes (noted in the supplement). 

We generated random (ER) graphs of size 10, 20, 50 and 100 nodes, creat- 
ing five different graphs for each number of nodes. The graphs were generated 
as in our convergence experiments, and our comparison to Monte Carlo test- 
ing are shown in Figure [5] where we demonstrate experimentally that our 
algorithm produces a lower bound. Our algorithm was run until the stan- 
dard deviation of fixation probabilities for all vertices was 2.5 x 10 -6 . The 
Monte Carlo simulations were each set at 10, 000 runs. 

8. Application: Competition Among Neural Axons 

In recent work, [26] created a model for synaptic competition based on 
death-birth updating under neutral drift. They noted that the model aligns 
well with their empirical observations. In the model, the graph represents a 
synaptic junction and the nodes represent sites in the junction. For every two 
adjacent sites in the synaptic junction, there is an undirected edge between 
the corresponding two nodes in the graph. Hence, in- and out- degrees of 
each node are the same. Initially, there are K different axon types located 
in the junction configured in a manner where all sites are initially occupied 
by one axon type. At each time step, an axon occupying one of the sites is 
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Figure 5: Mean-time-to- fixation comparison between algorithm and simulation. Note that 
the y-axis is a logarithmic scale. 



eliminated - making the site open. The selection of the axon for elimination 
(death) is with a uniform probability. Hence, there is no bias in this model. 
Following the elimination of an axon, an adjacent axon grows into the site. 
The adjacent axon is selected with a uniform probability of the eliminated 
axon's neighbors. Hence, based on the results of this paper, we can provide 
the following insights into synaptic competition. 

1. After t axons are eliminatedQthe probability of any site being occupied 
by an axon of a certain type can be calculated directly by Theorem |8j 
Even though there are K axon types, this theorem still applies as it 
only considers the probability of a node being a mutant (resp. a site 
being a one of the K axon types). 

2. Using point [T] above, we can determine the expected number of axons 
of a given type after t axons being eliminated. 



1 Note that the number of axons eliminated corresponds directly to the number of 
timesteps in the model. 
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3. After t axons are eliminated, the probability of any set of sites being 
occupied by a certain axon type is simply the sum of the probabilities 
of the individual sites being occupied by that axon. As a result, the 
fixation probability is additive. 

4. Leveraging point [3] above combined with an easy modification of the 
result of Broom et al. [6] for the BD model, the fixation probability 
of an axon originating at site i is ^ where ki is the number of sites 
adjacent to site i (hence the degree of node i in the corresponding 
graph) and is the total number of adjacencies in the synapse (hence, 
half the number of directed edges in the corresponding graph). 

5. Based on item [4] above and the results from Section [3j we can conclude 
that for a given axon type (let's call it "axon type A") occupying a 
set of sites, that if the average adjacencies of those sites is greater 
than (resp. less than) the overall average adjacencies for the sites in 
the entire synaptic junction, then as the number of eliminated axons 
approaches infinity, we can expect the number of axon type A in the 
synaptic junction will increase (resp. decrease) in expectation. 

6. We can directly apply Theorem [7] to find a lower-bound on the number 
of eliminated axons before fixation occurs. 

We note that the results stated above are either precise mathematical ar- 
guments or calculations that can be found exactly with a deterministic al- 
gorithm. They are not theoretical approximations and do not rely on sim- 
ulation. As such is the case, we can make more precise statements about 
synaptic competition (given the model) and can avoid the variance that 
accompanies simulation results. Insights such as these may lead to future 
biological experiments. 

9. Related Work 

Evolutionary graph theory was originally introduced in [TJ]. Previously, 
we have compiled a comprehensive review [21] for a general overview of the 
work in this exciting new area. 

While most work dealing with evolutionary graphs rely on Monte Carlo 
simulation, there are some good analytical approximations for the undi- 
rected/unweighted cased based on the degree of the vertices in question. 
Antal et al. [1] use the mean-field approach to create these approximations 
for the undirected/unweighted case. Broom et al. [B] derive an exact analyt- 
ical result for the undirected/unweighted case in neutral drift, which agrees 
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with the results of Antal et al. They also show that fixation probability is 
additive in that case (a result which we extend in this paper using a different 
proof technique). However, the results of [16] demonstrate that mean-field 
approximations break down in the case of weighted, directed graphs. [15] 
also studied weighted, directed graphs, but does so by using Monte Carlo 
simulation. [22] derive exact computation of fixation probability through 
means of linear programming. However, that approach requires an exponen- 
tial number of both constraints and variables and is intractable. The recent 
work of [2Z] introduces a parameter called graph determinacy which mea- 
sures the degree to which fixation or extinction is determined while starting 
from a randomly choses initial configuration. This property is then used 
to analyze some special cases of evolutionary graphs under birth-death up- 
dating. There has been some work on algorithms for fixation probability 
calculation that rely on a randomized approach [H [9] . [I] present a heuris- 
tic technique for speeding up Monte Carlo simulations by early termination 
while [5] present utilize simulation runs in a fully-polynomial randomized 
approximation scheme. However, our framework differs in that it does not 
rely on simulation at all and provides a deterministic result. Further, our 
non-randomized approach also allows for additional insights into the evolu- 
tionary process - such as monitoring the expected number of mutants as a 
function of time. Recently, [12] study the related problem of determining the 
probability of fixation given a single, randomly placed mutant in the graph 
where the vertices are "islands" and there are many individuals residing on 
each island in a well-mixed population. They use quasi-fixed points of ODE's 
to obtain an approximation of the fixation probability and performed experi- 
ments with a maximum of 5 islands (vertices) containing 50 individuals each. 
This continuous approximation provides the best results when the number 
of individuals in each island is much larger than the number of islands. As 
the problem of this paper can be thought of as a special case where each 
island has just one individual, it seems unlikely that the approximation of 
Houchmandzadeh and Vallade's approach will hold here. 

Some of the results in this paper were previously presented in conferences 
by the authors [23j E]- The analysis and experiments concerning the ex- 
pected number of mutants at a given time, the extension of the framework 
for other update rules (beyond birth-death), the use of the framework for the 
case of r > 1, and the neurology applications are all new results appearing 
for the first time in this paper. 
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10. Conclusion 

In this paper, we introduced a new approach to deal with problems relat- 
ing to evolutionary graphs that rely on "vertex probabilities." Our presented 
analytical method is the first deterministic method to compute fixation prob- 
ability and provides a number of novel uses and results for EGT problems: 

• Our method can be used to solve for the fixation probability under 
neutral drift orders of magnitude faster than Monte Carlo simulations, 
which is currently the presiding employed method in EGT studies. We 
have extended the method to all of the commonly used update pro- 
cesses in EGT. The special case of neutral drift is not only of interest 
in the literature [HI [15] but also it has been applied to problems in 
neurology [26] . 

• While the presented method is currently constrained to the case of 
neutral drift, we have demonstrated how it can inform cases of non- 
neutral drift by using it to provide both a lower and upper bound for 
this case. Combined with our analytical method's speed, this means 
that it can be used to acquire useful knowledge to guide general EGT 
studies interested in the case of advantageous mutants. 

• We have shown how our analytical method can be used to calculate a 
non-trivial lower bound to the mean time to fixation, providing a first 
step for a general method to computing this and related quantities that 
is lacking in the current literature. 

• We have shown how our method can be used to calculate deterministi- 
cally the expected number of mutants, which is useful for applications 
that require predictions on the number of mutants in the population 
under a specific finite time horizon. We have also provided results on 
the expected number of mutants on different common graph topologies, 
showing differences in the growth trajectories of amplifiers and sup- 
pressors on these different topologies. These results may prove highly 
significant in the recent application of EGT to distributed systems [13] 
where the problem of information diffusion is considered among com- 
puter systems. In such a domain, it may be insufficient to guarantee 
fixation in the limit of time - which may be impractical - but rather to 
make guarantees on the outcome of the process after a finite amount 
of time. 
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• Finally, we have shown how our method can provide insight when ap- 
plied to the problem of synaptic competition in neuroscience. 

Though evolutionary graph theory is still a relatively new research area, it 
is actively being studied in a variety of disciplines [141 fM\ 1281 1251 120] |26| 113] . 
We believe that more real-world applications will appear as this area gains 
more popularity As illustrated by recent work [221 HSj, experimental sci- 
entists with knowledge of EGT may be more likely to recognize situations 
where the model may be appropriate. As these cases arise, deterministic 
methods for addressing issues related to EGT may prove to be highly useful. 
However, this paper is only a starting point - there are still many impor- 
tant directions for future work. Foremost among such topics are scenarios 
where the topology of the graph also changes over time or where additional 
attributes of the nodes/edges in the graph affect the dynamics. 
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Supplementary Material 



11. Notes 

Throughout this supplement, we will use an extended notation. Fixation 
probability given initial configuration C is denoted Fc- For vertex % at time 
t, we denote this as Pr(M^). We will use Sf^ to denote the event that vertex 
i was selected for reproduction and to denote the event of i replacing 

j. We will often use conditional probabilities. For example, Pr(Mf^|C^) 
is the probability that Vi is a mutant given the initial set C of mutants. 
Throughout this supplement, unless noted otherwise, all of our probabili- 
ties will be conditioned on C^°\ We will drop it for ease of notation with 
the understanding that some set C of V were mutants at t — 0. Hence, 
Pr(Mf ) ) = Pr(Mf ) |C(°)). 

12. Proof of Theorem 1 



Vi, /zm^ 00 Pr(Mf ) |C(°)) = F c 
Proof. Consider the following definition property of Pr(M-^|C^ ^) 

Pr(Mf ) |C (0) ) = Pr(^ /W |C (0) ) (16) 

V'e2 v 
s.t. viev 

We note that as time approaches infinity, for all V E 2 V — — V we have 
Pr(V'M\CW) = 0. As Vi <£ 0, the statement follows. Q.E.D. 

13. Proof of Theroem 2 

Pr(M? } ) = 

Pr(M?- 1} ) + J2 w ^ • Pr(M5 t_1) ) ■ Pr(Sf |M? _1) ) - Wji ■ Pr(Mf _1) ) • Pr(Sf |M? _1) ) 

(vj,vi)eE 

Where Sf ^ is true iff Vi is selected for reproduction at time t. 
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Proof. Note we use the variable Rj i is true iff Vj replaces v% at time t. 
CLAIM 1: 

Pr(Mf ) ) = Pr(Mf" 1) A f\ -Sf) + £ Pr(Sf A R$ A M^) + 

J2 Pr(Sf A A Mf-)) 

(vj,vi)eE 

This is shown by a simple examination of exhaustive and mutually exclusive 
events based on the original model of 



CLAIM 2: 

Pr(Mf" 1) A /\ = Pr(Mf" 1) ). (l- £ Pr(Sf|M?-^ 

(Proof of claim 2) By exhaustive and mutual exclusive events, we have the 
following. 

Pr(Mf- 1} A f\ -Sf) = Pr(Mf- 1} )- £ Pr(Sf A M?^) 

(i)j,Di)e-E (vj,Vi)£E 

By the definition of conditional probability, we have the following 

Pr(M?- 1} A /\ = Pr(Mf- 1} )- £ (Pr(Sf IM^) ■ Pr(Mf- 1} )) 

= Pr(Mf" 1) ) -Pr(Mf _1) ) • Pr(Sf|Mf~ 1} ) 

(vj,Vi)£E 

= Pr(Mf^)(l- £ Pr(Sf| M r>; 

\ (vj,vi)eE 

The claim immediately follows. 

CLAIM 3: For all edges (vj,Vi), we have the following. 

Pr(Sf A R<f) A M$ t_1) ) = ^•Pr(Mj t - 1) )-Pr(S5 t) |M^ 1) ) 
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(Proof of claim 3) The following is a direct application of the definition of 
conditional probability. 

Pr(Sf ) A R$? A M? _1) ) = Pr(R| A M^ 1} |Sf ) • Pr(Sf) 

From our model, we note that given even the fitness of the nodes is not 

considered in determining if the event associated with is to occur. Hence, 

it follows that M^ _1) is independent of rJ? given Sf. As Pr(R| ) |Sj' ) ) = w jh 
we have the following. 

Pr(Sf A R^ A Mf- 1] ) = Pr(Rj t . ) \sf) • Pr(M^ 1} |Sf ) • Pr(Sf) 

= ^•Pr(M^ 1) |Sf)-Pr(Sf) 

By Bayes Theorem, and that the model causes ViPr(Sf ) ) > 0, we have the 
following. 



Pr(SfAR'?AM( M) ) = Wj ,-Pr(SflM^)- V J m ^ -Pr(Sf) 
j ■> j J j ■> Pr(S ) 

= ^•Pr(Mf- 1) )-Pr(S5 t) |M5 t - 1) ) 
The claim follows immediately. 

CLAIM 4: For all edges (vj,Vi), we have the following. 

Pr(Sf A-Rg AMf- 1 )) = (l-^)-Pr(Mf- 1) )-Pr(Sf|Mf- 1) ) 
(Proof of claim 4) This mirrors claim 3. 

(Proof of theorem) From claims 1-4, we have the following. 

Pr(Mf } ) = Pr(Mf (l- £ Pr^lM?" 1 ') J + (17) 



fa ■ Pr(Mj t - 1) ) ■ Pr(Sf|M( i - 1) )) + (18) 

(yj,vi)eE 

((l-w ji ).Pv(Mt 1) )--Pr(Sf\Mt 1) )) (19) 

{vj,vi)eE 

Which, after re-arranging some terms, gives us the statement of the theorem. 
Q.E.D. 
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14. Proof of Theorem 3 

When r = 1, if for some time t, Mi, the value Pr(M-^) is the same, then 
Pr(Mf ) ) = F c . 

Proof Sketch. Consider Pr{uf) = Pr(Mf _1) )+^ £(^)e£ w^iPriU^)- 
Pr(Mf _1) )) when fort -I, Mi,j we have Pr(M { -~ 1] ) = Pr(Mf~ 1} ). aear/y, 
in this case, the value for Pr(M^) = Pr(Mf -1 ^). ^4s the probabilities of all 
vertices was the same at t — 1, they remain so at t. Therefore, in this case, 
limt^Priuf) = Pr(Mf } ). QED 



15. Proof of Inequality 4 

For any time t, under neutral drift (r = 1), 

minPr(Mf } ) < F c < maxPr(Mf } ) 

i ' i 

Proof. PART 1: For any time t, under neutral drift (r = 1), Fc < 
maXiPr(Mf } ). 

We show that for each time step t, maxj Pr(M- t_1 - ) ) > maxj Pr(M^). Hence, 
by showing that, for any time t', we have maxj Pr(l\/lf ^ > lim t ^oo maxj Pr(M^) 
which by allows us to apply Theorem 1 and obtain the statement of this theo- 
rem. Suppose BWOC that at time t we have maxf Pr(M^~ 1 - ) ) < maxj Pr(M f^). 
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Then we have: 



maxPr(Mf 1) ) < 1 £ ' (Pr(M^) - Pr(Mf^)) 



(yj,vi)eE 
+Pr(Mf- 1} : 



< E ^•(maxPr(Mf 1) )-Pr(Mf- 1) )) 



(yj,Vi)eE 
,(t-ih 



+Pr(Mr 
g^f^ (mBxPr(Mr>) -Pr(Mr>)) 



maxPr(Mi ) < Pr (M^ LL -^ ) 

maxPr(Mf~ 1) ) < Pr(Mf~ 1} ) 

Which is clearly a contradiction and completes this part of the proof. 
PART 2: For any time t, under neutral drift (r = 1), Fq > minj Pr(M f 1 ). 
We show that for each time step t, minj Pr(Mf ^) < minj Pr(M^). Hence, 
by showing that, for any time t', we have minj Pr(M f < lirrit^^ maxj Pr(M^) 
which by allows us to apply Theorem 1 and obtain the statement of this theor- 
erm. Suppose BWOC that at time t we have mhi£ Pr(Mf 1} ) > minjPr(Mf ) ). 
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Then we have: 



minPr(Mf 1) ) > 1 £ ' (Pr(M^) - PrfM?- 1 ')) 



> 



(vj,Vi)£E 

+Pr(M? _1) ) 

1 £ ^,-(minPr(Mf 1) )-Pr(Mf- 1) )) 

(vj,Vi)eE 

+Pr(Mf~ 1} ) 



AT 

+Pr(Mf~ 1) 



- (^minPr(M^ 1) ) - Pr(Mf _1) )) 



minPr(Mf- 1) )(l- Efa '^ eg ^ ) > Pr(M^ } )(l - ^^^) 
minPr(M^ 1) ) > Pr(Mf _1) ) 

Which is clearly a contradiction and completes this part of the proof. Q.E.D. 

16. Proof of Theorem Theorem 4 

When r = 1 for disjoint sets C, D C V, F c + F D = F CuD . 
Proof. Consider some time t and vertex Vi. Clearly, by Corollary 1, Pr(IVlf ^) 
can be expressed as a linear combination of the form ^2 v . e y(Cj ■ Pr(M^)) 
where Cj is a coefficient. We note that these coefficients are the same re- 
gardless of the initial configuration of mutants that is conditioned on. 
Hence, Pr(Mf ) |C , (°)) is this positive function with Pr(Mj 0) ) = 1 if Vj G C 
and otherwise (see Proposition 3). Hence, for disjoint C, D, for any Vi G V, 
wehavePr(Mf ) |C(°)) + Pr(Mf ) | J D(°)) = Pr(Mf ] \(CUD)^). The statement 
follows. Q.E.D. 

17. Proof of Equation 9 



Pr(Mf')=(l-I).Pr(Mr)) + 1 £ Pr(M<-") 
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Under death-birth dynamics with neutral drift (r = 1). 

Proof. Df } and Bf } are random variables associated with birth and death 
events for vertex V{. 

CLAIM 1: Pr(M<*>) = Pi^M?^ A^Df } ) + £ ( ^ )6£ Pr^ AB « AM^) 
Follows directly form exhaustive and mutually exclusive events. 
CLAIM 2: Pr(Mf _1) A -Df } ) = (l - • Pr(Mf _1) ) 

By the definition of conditional probabilities, we have Pr(-iDf' ) |Mf _1 - ) ) • 
Pr(Mf 1 - > ). Also, we know the probability of a given node dying is always 
1/N. Hence, Pr(^Df } |Mf = Pr(^Df } ) = 1 - ^ and the claim follows. 
CLAIM 3: For any (vj,Vi) G E, we have 
Pr(D« A B« A Mf^) = • Pr(Mf x) ) 

in 

As both birth and death events occur independent of any node being a mutant 
at the previous time step, the definition of conditional probabilities gives us 
the following: 

Pr(D? ) A Bf A M$* _1) ) = Pr(D? } A Bf ) • Pr(M^ 1} ) (20) 

= Pr(Bj* ) |Df ) ) • Pr(Df ) ) • Pr(M^ 1) ) (21) 

From the model, we have the following: 

Pr(Bf|Df ) ) = (22) 
Pr(Df ) ) = 1/N (23) 

Hence, the claim follows. QED Q.E.D. 
18. Proof of Equation 10 



Pr(M? >) = ( 1 - *f ) • Pr(Mr>) + ^ £ Vv{M^) 



(vj,vi)eE 



\E\J v 1 ' \E i 
Under link dynamics with neutral drift (r = 1). 

Proof. Here S fj is the random variable associated with the selection of edge 

(Vi,Vj). 

CLAIM 1: Pr(Mf) = Pi^M?^ A f\ {vj , Vi)eE ^) + £ (iy>1li)6J5 Pr(S$ A 
M^) 
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Follows directly form exhaustive and mutually exclusive events. 
CLAIM 2: Pr(A (wjM ) 6B -S«) = l-f 

Clearly, we have Pr(A {v . ;V . )eE -Sj?) = Pr(V {(w)efWi} Sg). As there are 

fc£? incoming edges to ^, we know that Pr(V {( ^ )t)a)6£ |^ } S$) = 1 - ^, 
giving us the claim. 

CLAIM 3: Pr(Mf- 1) A A ( ^ )e£ -S<?) = (l - ff ) • P^M?^) 

For any a, /3, the random variable is independent from Mf 1 \ Hence, 
the claim immediately follows from this fact and claim 2. 
CLAIM 4: Pr(S| } A M$ t_1) ) = ± ■ Pr(M^ 1) ) 

As, by the definition of the model, Pr(Sj ) |Mj'" 1) ) = Pr(Sf?) = ±, the 
claim follows directly form the definition of conditional probabilities. QED 
Q.E.D. 

19. Proof of Theorem 5 

For a given set C, let F^\C) be the fixation probability under neutral 
drift and (C) be the fixation probability calculated using a mutant fitness 
r > 1. Then, under BD-B, BD-D, DB-B, DB-D, or LD dynamics, F«(C) < 
FW(C). 

Proof. First, some notation. 

• We define an interpretation, I : 2 y — > [0, 1] as probability distribution 
over mutant configurations. Hence, for some I we have ^2 v , e2 v 1(^0 = 
1. 

• Next, we define a transition function that maps configurations of mu- 
tants to probabilities, x '■ 2 y — > [0, 1] where for any C G 2 V , Xlc'e2 v C") — 
1. We will use x+ an d X- to indicate if the transition is made with a 
mutant being selected for birth (%+) or resident (%-). Hence, for some 

C G V and v i C, X-(C,CU{v}) = and X +(CU{v},C) = 0. Hence, 
for all C G 2 V , Z& & v(x+(C,C) + X -(C,C')) = 1. 

• If the transitioon function is based on birth-death and computed with 
some r > 1, then we will write it as x+\ X*- respectively. If computed 
with r = 1, then we write x+ d \ X'-^ respectively. 
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• For some C G 2 C , let inc(C) be the set of all elements D E 2 V s.t. 
|£>| > \C\ and x+{C,D) > 0. 

• For some C G 2 C , let dec(C) be the set of all elements D e 2 V s.t. 
\D\ < \C\ and X-(C,D) > 0. 

• Given set C C V, we will use to denot the probability of fixation 
given initial set of mutants C where the value r is used to calculate all 
transition probabilities. 

CLAIM 1: If a some time period, the probability distribution over mutant 
configurations is I, the fixation probability is ^2 Ce2 v 1(C) ' Fq\ 

Clearly, for any time t, = limi- > . 00 Pr(y^\C^). Under the assump- 
tion that there exists some tim u s.t. fixation is reached, we have: 

= Pr(V {uj) \C^) 
Pr(C(*)) 

Hence, F^ ] • Pr(C w ) = Pr(V (cj) A (?(*)). The statement then follows by the 
summation of exhaustive and mutually exclusive events. 

CLAIM 2: If a some time period t, the probability distribution over mutant 
configurations is I, and the transition functions used to reach the next time 
step are x+i X—> then the probability of being in some mutant configuration 
C at time t + 1 is given by £ D62 , 1(D) • C) + C)). 

Follows directly from the rules of dynamics. 

CLAIM 3: If a some time period t, the probability distribution over mu- 
tant configurations is I, mutant fitness r, and the transition functions used 

(r) (r) 

to reach the next time step are x+ > X- > an d a ^ subsequent transitions are 
computed using the same dynamics with neutral drift, then the fixation prob- 
ability is: 

V(I,r) = ^I(C). Yl (X ( ;\C,D).F^)+ (X { -\C,D)-FP)) 

Ce2 v \Deinc(C) D£dec(C) 
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Follows directly from claims 1-2. 



CLAIM 4: Under BD-B, BD-D, DB-B, DB-D, or LD dynamics, for some 

' < '< 

X ( -\C,D). 



r < r', for all C, D G 2 V , we have x+\C, D) < X + \C, D) and X-\C, D) > 



CLAIM 4a: For some r < r', for all C,D G 2 y , we have x+^C,^) < 

Let {vj} = D — C . For each vertex Vi, fi — 1 if i>j ^ C (a resident) and fi — r 
if i>j G C (a mutant). When D = C, the following are all summed over the 
set {fj G C|3i>i G C A (fi,^j) G £}. 

• Under BD-B, 

(r) (n n x \ ^ r • Wij 



tf<P-°) = E r .| C | + J V-| C j 



• Under BD-D, 



N ■ J2 Vq \(vi,v q )eE w ig ' fg 1 



Vi eC\ t-'Vq^i^q 

• Under DB-B, 

(vi,vj)eE 

• Under DB-D, 

XJ'(C,D) = V 

(vi,vj)eE 

• Under LD, 

x<?(c,d) = y: v — — — r 

vieCKv^eE ^v q ,v e \(.v q ,v e )eE w ge ' h 
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By simple algebraic manipulation, for each of these, when all values other 
than r are fixed, they increase as r increases. 



CLAIM 4b: For some r < r' , for all C, D £ 2 V , we have X - (C,D) > 
X-\C,D). Let {vj} = C — D. For each vertex v i: fa — 1 if Vi ^ C (a 
resident) and /, = r if vi £ C (a mutant). When D = C, the following are 
all summed over the set {vj £ V — C\3v,i £ V — C A (v i: Vj) £ E}. 

• Under BD-B, 



x { ;\c,d) 



E 



r-\C\+N-\C 



v l ev-c\ 
(vi,vj)eE 



• Under BD-D, 



ViEV-C\ 

(vi,vj)eE 



E 




• Under DB-B, 



X ( l\C,D) 



E 



iv-E. 



(vi,vj)eE 



• Under DB-D, 



xV r, (c,^) 



ViEV-C\ 

(vi,vj)eE 



E 




• Under LD, 



^6y-c|(^,^)eB 



E 
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By simple algebraic manipulation, for each of these, when all values other 
than r are fixed, they decrease as r increases. 

CLAIM 5: Given some C G 2 V , for all pairs D, D' where D G inc(C) and 
jy G dec(C), we have > F§) . 

Follows directly from Theorem 5. 

CLAIM 6: Given interpretation I, under BD-B, BD-D, DB-B, DB-D, or LD 
dynamics, for some r > 1, V(I, r) > V(I, 0). 

Let us consider some set C from the outermsot summation in the computa- 
tion of V(I, r). Suppose, BWOC, there exists some C G 2 V s.t.: 

J2 (XV(C,D).F$))+ £ (X W (C,D).F«) < J2 (x't ) (C,D).F^) + £ ( X ^(C,D)-F^) 

Deinc(C) Dedec(C) Deinc(C) Dedec(C) 

This give us: 

£ (x«(C,D).FW)- £ ( X W(C7,D).FW) < £ ( X W(C,D). F W)- J2 (xL r) (C,D).F«) 

D6int(C) D£inc(C) Dedec(C) Dedec(C) 

Let F sm = ira/XF^I-D G mc(C)} and = sup{F^ ] \D G dec(C)}, this give 
us: 

F 3m £ (x^ ) (C,D)- x ^ 1) (C,D)) < F i9 £ (xL 1) (C,D)- x L r) (C,D)) 

D£mc(C) Dedec(C) 

Consider ther following: 

J2 X { ;\C,D)+ J2 X M (C,D) = J2 X { 1\C,D)+ Y. X W (C,D) 

DGinc(C) Dedec(C) Deinc(C) Dedec(C) 

Deinc(C) Deinc(C) Dedec(C) Dedec(C) 

Note that by claim 4, both sides of the above equation are positive numbers. 
Hence, we have F sm < Fi g , which contradicts claim 5. 

PROOF OF THEOREM: Let P (1) (I,r) = P(I,r) and V {i+1) (I, r) = V(V {i) (I,r),r). 
By claim 6, for any i, V^ l+1 \l,r) > V(V^(I,r),r). Consider an inter- 
pretation I that describes the initial probability distribution over mutant 
configurations. The fixation probability under neutral drift is V(I,1). For 
some value r > 1, the fixation probaiblity is £wn i _ KX ,'P^(I, r). Clearly, 
Um ir + O0 V®{l,r)>V{l,l)- QE.D. 
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20. Proof of Theorem 6 



tc = jrJ2 t -^'- F ^ 

c t=i 



Proof. This proof was first presented in [2]. The mean time to fixation is 
described as the expected time to fixation given that the process fixates. Let 
F t \c be the probability that fixation is reached in exactly t time-steps or less. 
Hence, the probability of reaching fixation in exactly t time steps conditioned 
on the process reaching fixation is {F t \c — Ft-\\c)/Fc- The remainder of the 
theorem follows from the definition of an expected value. Q.E.D. 

21. Proof of Theorem 7 

We introduce two pieces of notation tfi X ,t convg . We define tf ix as a time 

s-t. t c = ^E&* • {Ft\c ~ Ft-i\c). and t convg s.t. Vi, Pr(Mf— } ) = F c . 
While in reality, both of these values could be infinite, we note that the 
relationship oo > tfi X > t convg holds and that using a lower value for tf ix 
and/or t convg will still ensure we have a lower bound. 



j oo ^ oo 

JT/J ' {Prwn,t - •Pmin.t-l) < JT t ' ( F t\C ~ F t-l\C') 
° t=l C t=l 

Where P minj4 = min i Pr(Mf ) ). 
Proof. First, we have the following. 

tconvg tfix 

1 ■ ( P ™*J - P min,t-l) < £ t- (F t]C - F t _ 1]C ) (24) 

t=l t=l 

For any time t, let V { § = F t[c - F t _ 1]c and Pr(AM^ in ) = ^U,t - Pmm,t-i- 
As for each time t, we have P m m,t > Ft\c, we can define 0® as the portion of 
Vq accounted for in Pr(AM^). This results in having 9# = whenever 
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t' > t or t > t convg as well as the following: 



1 c — 



(t) s 
min/ 



Pr(AM 

tconvg 

^ ^ t ' (yPmm,t -Pmin,t— 1 



t=l 



£«<;> 

t'=i 

t fix 

t'=t 

^convg ^fix 

= EE« 

t=i t'=t 

tconvg tfix 

= EE'« 

t=i t'=i 

tfix tconvg 

= EE« 

t'=i t=i 

tfix t! 

= EE« 

t'=i t=i 



(*') 



(*') 



(*') 



(f) 



We also note that the following is true: 

tfix 



y^j ■ ( f ac - F t -i\ C ) 



t=l 



tfix 

t'=i 

tfix t! 



= E>'E> 



(*') 



t'=i t=i 
t fix t> 



>- J2J2 te 



it') 



t'=i t=i 

tconva 



t ■ (Pmin,t ~ ^min,t-l) 



t=l 



Which concludes the proof. Q.E.D. 



22. Materials and Methods 



(25) 
(26) 
(27) 
(28) 
(29) 
(30) 



(31) 
(32) 
(33) 
(34) 



Except for the experiments dealing with time to fixation/extinction, all 
algorithms were implemented in Python and run on a 2.33GHz Intel Xeon 
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CPU. The time-to-fixation experiments were run on a machine equipped with 
an Intel Core i7 M620 processor running at 2.67 GHz with 4 GB RAM. 

All graphs in the experiments were generated using the Python NetworkX 
package [H agberg et al. . Parameters used for the experiments concerning the 



expected number of mutants were m=l for BA, p = 0.5 for ER, and k = 2 
and p = 0.5 for NWS graph generator functions. 

We modified Algorithm 1-ACC based on the results on mean time to 
fixation as follows: 

• Before line 14, insert: t += 1; Sum += t*(min(p)-min(q)) 

• Replace line 17 with: return = Sum/average(p), where average(q) 
is the algorithm's best estimate for P c at termination. 
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